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Abstract 

A traditional assumption underlying most data converters is that the signal should be sampled at a rate exceeding 
twice the highest frequency. This statement is based on a worst-case scenario in which the signal occupies the 
entire available bandwidth. In practice, many signals are sparse so that only part of the bandwidth is used. In this 
paper, we develop methods for low-rate sampling of continuous-time sparse signals in shift-invariant (SI) spaces, 
generated by m kernels with period T. We model sparsity by treating the case in which only k out of the m 
generators are active, however, we do not know which k are chosen. We show how to sample such signals at a rate 
much lower than mAT, which is the minimal sampling rate without exploiting sparsity. Our approach combines ideas 
from analog sampling in a subspace with a recently developed block diagram that converts an infinite set of sparse 
equations to a finite counterpart. Using these two components we formulate our problem within the framework of 
finite compressed sensing (CS) and then rely on algorithms developed in that context. The distinguishing feature 
of our results is that in contrast to standard CS, which treats finite-length vectors, we consider sampling of analog 
signals for which no underlying finite-dimensional model exists. The proposed framework allows to extend much 
of the recent literature on CS to the analog domain. 

I. Introduction 

Digital applications have developed rapidly over the last few decades. Signal processing in the discrete domain 
inherently relies on sampling a continuous-time signal to obtain a discrete-time representation. The traditional 
assumption underlying most analog-to-digital converters is that the samples must be acquired at the Shannon- 
Nyquist rate, corresponding to twice the highest frequency [1], [2]. 

Although the bandlimited assumption is often approximately met, many signals can be more adequately modeled 
in alternative bases other than the Fourier basis [3], [4], or possess further structure in the Fourier domain. Research 
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in sampling theory over the past two decades has substantially enlarged the class of sampling problems that can 
be treated efficiently and reliably. This resulted in many new sampling theories which accommodate more general 
signal sets as well as various linear and nonlinear distortions [5], [6], [4], [7], [8], [9], [10], [11]. 

A signal class that plays an important role in sampling theory are signals in shift-invariant (SI) spaces. Such 
functions can be expressed as Unear combinations of shifts of a set of generators with period T [12], [13], [14], 
[15], [16]. This model encompasses many signals used in communication and signal processing. For example, 
the set of bandlimited functions is SI with a single generator. Other examples include splines [4], [17] and pulse 
amplitude modulation in communications. Using multiple generators, a larger set of signals can be described such 
as multiband functions [18], [19], [20], [21], [22], [23]. Sampling theories similar to the Shannon theorem can be 
developed for this signal class, which allows to sample and reconstruct such functions using a broad variety of 
filters. 

Any signal x{t) in a SI space generated by m functions shifted with period T can be perfectly recovered from 
m sampling sequences, obtained by filtering x{t) with a bank of m filters and uniformly sampling their outputs 
at times nT. The overall sampling rate of such a system is m/T. In Section Ull we show explicitly how to recover 
x{t) from these samples by an appropriate filter bank. If the signal is generated by k out of the m generators, then 
as long as the chosen subset is known, it suffices to sample at a rate of k/T corresponding to uniform samples 
with period T at the output of k filters. However, a more difficult question is whether the rate can be reduced if 
we know that only k of the generators are active, but we do not know in advance which ones. Since in principle 
x{t) may be comprised of any of the generators, it may seem at first that the rate cannot be lower than m/T. 

This question is a special case of sampling a signal in a union of subspaces [24], [25], [26]. In our problem, 
x{t) lies in one of the subspaces expressed by k generators, however we do not know which subspace is chosen. 
Necessary and sufficient conditions where derived in [24], [25] to ensure that a sampling operator over such a 
union is invertible. In our setting this reduces to the requirement that the sampling rate is at least 2k /T. However, 
no concrete sampling methods where given that ensure efficient and stable recovery, and no recovery algorithm was 
provided from a given set of samples. Finite-dimensional unions where treated in [26], for which stable recovery 
methods where developed. Another special case of sampling on a union of spaces that has been studied extensively 
is the problem underlying the field of compressed sensing (CS). In this setting, the goal is to recover a length m 
vector X from p < m linear measurements, where x is known to be /c-sparse in some basis [27], [28]. Many stable 
and efficient recovery algorithms have been proposed to recover x in this setting [29], [30], [31], [32], [28], [33], 
[26]. 

A fundamental difference between our problem and mainstream CS papers is that we aim to sample and 



reconstruct continuous signals, while CS focuses on recovery of finite vectors. The methods developed in the 
context of CS rely on the finite nature of the problem and cannot be immediately adopted to infinite-dimensional 
settings without discretization or heuristics. Our goal is to directly reduce the analog sampling rate, without first 
requiring the Nyquist-rate samples and then applying finite-dimensional CS techniques. 

Several attempts to extend CS ideas to the analog domain were developed in a set of conferences papers [34], 
[35]. However, in both papers an underlying discrete model was assumed which enabled immediate application of 
known CS techniques. An alternative analog framework is the work on finite rate of innovation [36], [37], in which 
x{t) is modeled as a finite linear combination of shifted diracs (some extensions are given to other generators 
as well). The algorithms developed in this context exploit the similarity between the given problem and spectral 
estimation, and again rely on finite dimensional methods. 

In contrast, the model we treat in this paper is inherently infinite dimensional as it involves an infinite sequence 
of samples from which we would like to recover an analog signal with infinitely many parameters. In such a setting 
the measurement matrix in standard CS is replaced by a more general linear operator. It is therefore no longer clear 
how to choose such an operator to ensure stability. Furthermore, even if a stable operator can be implemented, 
it will result in infinitely many compressed samples. As standard CS algorithms operate on finite-dimensional 
optimization problems, they caimot be applied to infinite dimensional sequences. 

In our previous work, we considered a sparse analog samphng problem in which the signal x{t) has a multiband 
structure, so that its Fourier transform consists of at most N bands, each of width limited by B [21], [22], [23]. 
Explicit sub-Nyquist sampling and reconstruction schemes were developed in [21], [22], [23] that ensure perfect 
recovery of multiband signals at the minimal possible rate, without requiring knowledge of the band locations. 
The proposed algorithms rely on a set of operations grouped under a block named continuous-to-finite (CTF). 
The CTF, which is further developed in [38], essentially transforms the continuous reconstruction problem into a 
finite dimensional equivalent, without discretization or heuristics. The resulting problem is formulated within the 
framework of CS, and thus can be solved efficiently using known tractable algorithms. 

The samphng methods used in [21], [23] for bhnd multiband sampling are tailored to that specific setting, and 
are not applicable to the more general model we consider here. Our goal in this paper is to capitalize on the key 
elements from [21], [38], [23] that enable CS of multiband signals and extend them to the more general SI setting 
by combining results from standard sampling theory and CS. Although the ideas we present are rooted in our 
previous work, their apphcation to more general analog CS is not immediately obvious. To extend our work, it 
is crucial to setup the more general problem treated here in a particular way. Therefore, a large part of the paper 
is focused on the problem setup, and reformulation of previously derived results. We then show explicitly how 



signals in a SI union created by m generators with period T, can be sampled and stably recovered at a rate much 
lower than m/T. Specifically, if k out of m generators are active, then it is sufficient to use 2k < p < m uniform 
sequences at rate 1/T, where p is determined by the requirements of standard CS. 

The paper is organized as follows. In Section |ll] we provide background material on Nyquist-rate sampling in SI 
spaces. Although most of these results are known in the literature, we review them here since our interpretation of 
the recovery method is essential in treating the sparse setting. The sparse SI model is presented in Section |llll In 
this section we also review the main elements of CS needed for the development of our algorithm, and elaborate 
more on the essential difficulty in extending them to the analog setting. The difference between sampling in general 
SI spaces and our previous work [21] is also highlighted. In Section IVl we present our strategy for CS of SI analog 
signals. Some examples of our framework are discussed in Section |Vll 

II. Background: Sampling in SI Spaces 

Traditional sampling theory deals with the problem of recovering an unknown function x{t) G L2 from its 
uniform samples at times t = nT where T is the sampling period. More generally, the signal may be pre-filtered 
prior to sampling with a filter s*{—t) [4], [7], [11], [39], [16], where (•)* denotes the complex conjugate, as 
illustrated in the left-hand side of Fig. [T] The samples c[n] can be represented as the inner products c[n] = 



xit) 



s%-t) 



c[n\ 




d[n] 









t = nT 





a{t) 




) — ► 





xit) 



Fig. I. Non-ideal sampling and reconstruction. 
{s{t — nT),x{t)) where 
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In order to recover x{t) from these samples it is typically assumed that x{t) lies in an appropriate subspace A of 

A common choice of subspace is a SI subspace generated by a single generator a{t). Any x{t) G A has the 
form 

x(t) = ^d[n]a(t-nT), (2) 



for some generator a{t) and sampling period T. Note that d[n] are not necessarily pointwise samples of the signal. 
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where we defined 

then x(t) can be perfectly reconstructed from the samples c[n] in Fig. [T] [39], [6]. The function 4>sA{e^'^) is the 
discrete-time Fourier transform (DTFT) of the sampled cross-correlation sequence: 

rsa[n] = {s{t-nT),a{t)). (5) 

To emphasize the fact that the DTFT is 27r-periodic we use the notation ^{e^'^). Recovery is obtained by filtering 
the samples c[n] by a discrete-time filter with frequency response 

followed by modulation by an impulse train with period T and filtering with an analog filter a{t). The overall 
sampling and reconstruction scheme is illustrated in Fig. [U Evidently, SI subspaces allow to retain the basic flavor 
of the Shannon sampling theorem in which sampling and recovery are implemented by filtering operations. 

In this paper we consider more general SI spaces, generated by m functions a£{t) G L2,l < £ < m. A 
finitely-generated SI subspace in L2 is defined as [12], [13], [15] 

A = I x{t) =f2Yl di[n]ai{t - riT) : de[n] G ^2 } • (7) 

I £=1 n& ) 

The functions ae{t) are referred to as the generators of A. In the Fourier domain, we can represent any x{t) ^ A 



as 



= j;Z),(e^-^)A,M, (8) 

where 

^Ke^'") = E^^Ne^'"", (9) 

neZ 

is the DTFT of d^[n]. Throughout the paper, we use upper-case letters to denote Fourier transforms: X{uj) is the 
continuous-time Fourier transform of the function x{t), and C{e^^) is the DTFT of the sequence c[n\. 

In order to guarantee a unique stable representation of any signal in A by coefficients di[n\, the generators a£{t) 
are typically chosen to form a Riesz basis for L2. This means that there exists constants q > and /? < 00 such 



that 



"^"^ di[n\ai{t - nT) 
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where ||d|p = X^^^ Ylnez M^NP' ™d the norm in the middle term is the standard L2 norm. By taking Fourier 
transforms in ([TOb it follows that the generators a^(t) form a Riesz basis if and only if [13] 



(11) 



where 



MAAien 



(12) 



Here (pA^Aeie^'^) is defined by (Hi with Ai{e^'^), Ai{e^'^) replacing 5(e-''^), ^(e-?'^). Throughout the paper we 
assume that ([TT]) is satisfied. 

Since x{t) lies in a space generated by m functions, it makes sense to sample it with m filters as in the 

left-hand side of Fig. |2] The samples are given by 
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Fig. 2. Sampling and reconstruction in shift-invariant spaces. 



ci[n] = {s^{t — nT)^x{t)) , 1 < £ < m. 



(13) 



The following proposition provides a simple Fourier-domain relationship between the samples ci[n] of ( 1131 ) and 



the expansion coefficients de[n\ of (O: 

Proposition 1: Let Ci[n\ = — nT),x{t)), 1 < ^ < m be a set of m sequences obtained by filtering x{t) 
of dV) with m filters s}{—t) and sampling their outputs at times nT, as depicted in the left-hand side of Fig. |2] 
Denote by c{e^'^),d{e^'^) the vectors with ^th elements Ci{e^^), Di{e^'^), respectively. Then 



where 



c(e^''^) = M5A(e^''')d(e 
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and (^s^aX^^'^) ^re defined by (HJl. 

Proof: The proof follows immediately by taking the Fourier transform of ([T3l) : 
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where we used ([8]l and the fact that D{e^^) is 27r-periodic. In vector from, ([T6l ) reduces to ([141 ). ■ 
Proposition [T] can be used to recover x{t) from the given samples as long as 'WisA{e.^'^) is invertible a.e. 
in uj. Under this condition, the expansion coefficients can be computed as d(eJ'^) = M.'g\{ei'^)c{e^'^). Given 
di[n], x{t) is formed by modulating each coefficient sequence by a periodic impulse train X^ngz ^(^ ~ '^-^) ^^^^ 
period T, and filtering with the corresponding analog filter a^(t). In order to ensure stable recovery we require 
al ^ NLsAie-''^) ^ /3I a.e. In particular, we may choose Si{t) = ai{t) due to ([TT] ). The resulting sampling and 
reconstruction scheme is depicted in Fig. |2] 

The approach of Fig. |2] results in m sequences of samples, each at rate 1/T, leading to an average sampling 
rate of m/T. Note that from d?) it follows that in each time step T, x{t) contains m new parameters, so that the 
signal has m degrees of freedom over every interval of length T. Therefore this sampling strategy has the intuitive 
property that it requires one sample for each degree of freedom. 



III. Union of Shift-Invariant Subspaces 

Evidently, when subspace information is available, perfect reconstruction from linear samples is often achievable. 
Furthermore, recovery is possible using a simple filter bank. A more interesting scenario is when x{t) lies in a 



union of SI subspaces of the form [26] 
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\l\=k 

where the notation \£\ = k means a union (or sum) over at most k elements. Here we consider the case in which 
the union is over k out of m possible subspaces {Ae, I < i < m}, where is generated by a^(t). Thus, 



so that only k out of the m sequences ddn] in the sum ([TS] ) are not identically zero. Note that ([TS] ) no longer 
defines a subspace. Indeed, while the union contains signals of the form di[n\ai{t — nT) for two distinct values 
of i, it does not include their linear combinations. 

In principle, if we know which k sequences are non-zero, then x{t) can be recovered from samples at the output 
of k filters using the scheme of Fig. |2] The resulting average sampling rate is k/T since we have k sequences, 
each at rate 1/T. Alternatively, even without knowledge of the active subspaces, we can recover x{t) from samples 
at the output of m filters resulting in a sampling rate of m/T. Although this strategy does not require knowledge 
of the active subspaces, the price is an increase in sampling rate. 

In [24], [25], the authors developed necessary and sufficient conditions for a sampling operator to be invertible 
over a union of subspaces. Specializing the results to our problem implies that a minimal rate of at least 2k /T is 
needed in order to ensure that there is a unique SI signal consistent with the samples. Thus, the fact that we do not 
know the exact subspace leads to an increase of at least a factor two in the minimal rate. However, no concrete 
methods were provided to reconstruct the original signal from its samples. Furthermore, although conditions for 
invertibility were provided, these do not necessarily imply that a stable and efficient recovery is possible at the 
minimal rate. 

Our goal is to develop algorithms for recovering x{t) from a set of 2/;; < p < m sampling sequences, obtained 
by sampling the outputs of p filters at rate 1/T. Before developing our sampling scheme, we first explain the 
difficulty in addressing this problem and its relation to prior work. 

A. Compressed Sensing 

A special case of a union of subspaces that has been treated extensively is CS of finite vectors [27], [28]. In 
this setup, the problem is to recover a finite-dimensional vector x of length m from p < m linear measurements 
y where 




(18) 
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y = Mx, 



(19) 



for some matrix M of size p x m. Since the equations ([T9l ) are underdetermined, more information is needed in 
order to recover x. The prior assumed in the CS literature is that x = $a where $ is an m x m invertible matrix, 
and a is fc-sparse, so that it has at most k non-zero elements. This prior can be viewed as a union of subspaces 
where each subspace is spanned by k columns of 

A sufficient condition for the uniqueness of a fc-sparse solution a to the equations y = M$a is that A = M$ 
has a Kruskal-rank of at least 2k [32], [40]. The Kruskal-rank is the maximal number q such that every set of q 
columns of A is linearly independent [41]. This unique a can be recovered by solving the optimization problem 
[27]: 

niin||Q;||o s. t. y = Act, (20) 

where the pseudo-norm Iq counts the number of non-zero entries. Therefore, if we are not concerned with stability 
and computational complexity, then 2k measurements are enough to recover a exactly. Since (l20b is known to 
be NP-hard [27], [28], several alternative algorithms have been proposed in the literature that have polynomial 
complexity. Two prominent approaches are to replace the norm by the convex £i norm, and the orthogonal 
matching pursuit algorithm [27], [28]. For a given sparsity level k, these techniques are guaranteed to recover 
the true sparse vector a as long as certain conditions on A are satisfied, such as the restricted isometry property 
[28], [42]. The efficient methods proposed to recover x all require a number of measurements p that is larger than 
2k, however still considerably smaller than m. For example, if A is chosen as p random rows from the Fourier 
transform matrix, then the £i program will recover a with overwhelming probability as long as p > ck log m where 
c is a constant. Other choices for A are random matrices consisting of Gaussian or Bernoulli random variables 
[33], [43]. In these cases, on the order of k\og{m/k) measurements are necessary in order to be able to recover 
Q efficiently with high probability. 

These results have also been generalized to the multiple-measurement vector (MMV) model in which the problem 
is to recover a matrix X from matrix measurements Y = AX, where X has at most k non-zero rows. Here again, 
if the Kruskal-rank of A is at least 2k, then there is a unique X consistent with Y. This unique solution can be 
obtained by solving the combinatorial problem 

min|/(X)|, s.t. Y = AX, (21) 

where /(X) is the set of indices corresponding to the non-zero rows of X [40]. Various efficient algorithms that 
coincide with (12111 under certain conditions on A have also been proposed for this problem [44], [40], [3 8]. 



B. Compressed Sensing of Analog Signals 

Our problem is similar in spirit to finite CS: we would like to sense a sparse signal using fewer measurements 
than required without the sparsity assumption. However, the fundamental difference between the two stems from 
the fact that our problem is defined over an infinite-dimensional space of continuous functions. As we now show, 
trying to represent it in the same form as CS by replacing the finite matrices by appropriate operators, raises 
several difficulties that precludes direct application of CS-type results. 

To see this, suppose we represent x{t) in terms of a sparse expansion, by defining an infinite-dimensional 
operator <I>(t) corresponding to the concatenation of the functions a^(t — riT), and an infinite sequence a G £2 
which consists of the concatenation of the sequences di[n\. We may then write x{t) = ^{t)a which resembles 
the finite expansion x = Since dg[n\ is identically zero for several values of £, a will contain many zero 
elements. Next, we can define a measurement operator M (t) so that the measurements are given by y = Aa where 
A = M{t)<^>{t). 

In analogy to the finite setting, the recovery properties of a should depend on A. However, immediate application 
of CS ideas to this operator equation is impossible. As we have seen, the ability to recover ct in the finite setting 
depends on its sparsity. In our case, the sparsity of a is always infinite. Furthermore, a practical way to ensure 
stable recovery with high probability in conventional CS is to draw the elements of A at random, with the number 
of rows of A proportional to the sparsity. In the operator setting, we cannot clearly define the dimensions of A 
or draw its elements at random. Even if we can develop conditions on A such that the measurement sequence 
y = Aa uniquely determines a, it is still not clear how to recover a from y. The immediate extension of basis 
pursuit to this context would be: 

min||a||i s.t. y = Aa. (22) 

a 

Although ((22)) is a convex problem, it is defined over infinitely many variables, with infinitely many constraints. 
Convex programming techniques such as semi-infinite programming and generalized semi-infinite programming, 
allow only for infinite constraints while the optimization variable must be finite. Therefore, (l22l) cannot be solved 
using standard optimization tools as in finite-dimensional CS. 

This discussion raises three important questions we need to address in order to adapt CS results to the analog 
setting: 

1) How do we choose an analog sampling operator? 

2) Can we introduce structure into the sampling operator and still preserve stability? 

3) How do we solve the resulting infinite-dimensional recovery problem? 



C. Previous Work on Analog Compressed Sensing 

In our previous work [21], [22], [23], we treated a special case of analog CS. Specifically, we considered blind 
multiband sampling in which the goal is to sample a bandlimited signal x{t) whose frequency response consists 
of at most bands of length limited by B, with unknown support. In Section |Vl] we show that this model can 
be formulated as a union of SI subspaces ([TSl l. In order to sample and recover x{t) at rates much lower than 
Nyquist, we proposed two types of sampling methods: In [21] we considered multicoset sampling while in [23] 
we modulated the signal by a periodic function, prior to standard low-pass sampling. Both strategies lead to p 
sequences of low -rate uniform samples which, in the Fourier domain, can be related to the unknown x{t) via an 
infinite measurement vector (IMV) model [38], [21]. This model is an extension of the MMV problem to the case 
in which the goal is to recover infinitely many unknown vectors that share a joint sparsity pattern. Using the IMV 
techniques developed in [21], [38], x{t) can then be recovered by solving a finite-dimensional convex optimization 
problem. We elaborate more on the IMV model below, as it is a key ingredient in our proposed sampling strategy. 

The sampling method described above is tailored to the multiband model, and exploits the fact that the spectrum 
has many intervals that are identically zero. Applying this approach to the general SI setting will not lead to perfect 
recovery. In order to extend our previous results, we therefore need to reveal the key ideas that allow CS of analog 
signals, rather than analyzing a specific set of sampling equations as in [21], [22], [23]. Further study of blind 
multiband sampling suggests two key elements enabling analog CS: 

1) Fourier domain analysis of the sequences of samples. 

2) Choosing the sampling functions such that we obtain an IMV model. 

Our approach is to capitalize on these two components and extend them to the model ([TSl l. 

To develop an analog CS system, we design p < m sampling filters Si{t) that enable perfect recovery of x{t). In 
view of our previous discussion, our task can be rephrased as determining sampling filters such that in the Fourier 
domain, the resulting samples can be described by an IMV system. In the next section we review the key elements 
of the IMV problem. We then show how to appropriately choose the sampling filters for the general model ([TSl l. 

IV. Infinite Measurement Model 
In the IMV model the goal is to recover a set of unknown vectors x(A) from measurement vectors 

y(A) = Ax(A), A G A, (23) 

where A is a set whose cardinality can be infinite. In particular, A may be uncountable, such as the frequencies 
u) G [— 7r,7r). The /c-sparse IMV model assumes that the vectors {x(A)}, which we denote for brevity by x(A), 



share a joint sparsity pattern, that is, the non-zero elements are supported on a fixed location set of size k [38]. 

As in the finite-case it is easy to see that if (t(A) > 2k, where cr(A) is the Kruskal-rank of A, then x(A) is the 
unique A;-sparse solution of (l23l ) [38]. The major difficulty with the IMV model is how to recover the solution set 
x(A) from the infinitely many equations ( |23l ). One suboptimal strategy is to convert the problem into an MMV by 
solving (|23] ) only over a finite set of values A. However, clearly this strategy cannot guarantee perfect recovery. 
Instead, the approach in [38] is to recover x(A) in two steps. First, we find the support set S of x(A), and then 
reconstruct x(A) from the data y(A) and knowledge of S. 

Once S is found, the second step is straightforward. To see this, note that using S, (1231 ) can be written as 

y(A) = A5X^(A), AG A, (24) 

where A5 denotes the matrix containing the columns of A whose indices belong to S, and x'^(A) is the vector 
consisting of entries of x(A) in locations S. Since x(A) is A;-sparse, |5| < k. Therefore, the columns of As 
are linearly independent (because o"(A) > 2k), implying that aI^As = I, where A^ = (A^Ag) ^ A^ is the 
pseudo-inverse of A5 and (•)^ denotes the Hermitian conjugate. Multiplying (l24l ) by A^ on the left gives 

x^(A) = A^y(A), A G A. (25) 

The components in x(A) not supported on S are all zero. Therefore ( |25] ) allows for exact recovery of x(A) once 
the finite set S is correctly identified. 

It remains to determine S efficiently. In [38] it was shown that S can be found exactly by solving a finite MMV. 
The steps used to formulate this MMV are grouped under a block referred to as the continuous-to-finite (CTF) 
block. The essential idea is that every finite collection of vectors spanning the subspace span(y(A)) contains 
sufficient information to recover S, as incorporated in the following theorem [38]: 

Theorem 1: Suppose that o"(A) > 2k, and let V be a matrix with column span equal to span(y(A)). Then, the 
linear system 

V = AU (26) 
has a unique A;-sparse solution U whose support is equal 5. 

The advantage of Theorem [T] is that it allows to avoid the infinite structure of (1231 ) and instead find the finite set S 
by solving the single MMV system of (l26l) . The additional requirement of Theorem [T] is to construct a matrix V 
having column span equal to span(y(A)). The following proposition, proven in [38], suggests such a procedure. 
To this end, we assume that x(A) is piecewise continuous in A. 
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Fig. 3. The fundamental stages for the recovery of the non-zero location set 5* in an IMV model using only one finite-dimensional problem. 



Proposition 2: If the integral 
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y(A)y^(A)dA, 



(27) 



AeA 



exists, then every matrix V satisfying Q = VV^ has a column span equal to span(y(A)). 

Fig. [3j taken from [38], summarizes the reduction steps that follow from Theorem [T] and Proposition |2l Note, 
that each block in the figure can be replaced by another set of operations having an equivalent functionality. In 
particular, the computation of the matrix Q of Proposition |2] can be avoided if alternative methods are employed for 
the construction of a frame V for span(y(A)). In the figure, / indicates the joint support set of the corresponding 
vectors. 



V. Compressed Sensing of SI Signals 

We now combine the ideas of Sections |ll] and |IV] in order to develop efficient sampling strategies for a union 
of subspaces of the form (fTSl ). Our approach consists of filtering x{t) with p < m filters Si{t), and uniformly 
sampling their outputs at rate 1/T. The design of Si{t),l < i < p relies on two ingredients: 

1) A matrix A chosen such that it solves a discrete CS problem in the dimensions m (vector length) and k 
(sparsity). 

2) A set of functions h.i{t), \ < i < m which can be used to sample and reconstruct the entire set of generators 
(ii{t),l < i < m, namely such that 'M.uAi^^^'^) is stably invertible a.e. 

The matrix A is determined by considering a finite-dimensional CS problem where we would like to recover a 
fe-sparse vector x of length m from p measurements y = Ax. The value of p can be chosen to guarantee exact 
recovery with combinatorial optimization, in which case p > 2k, or to lead to efficient recovery (possibly only 
with high probability) requiring p > 2k. We show below that the same A chosen for this discrete problem can 
be used for analog CS. The functions hi{t) are chosen so that they can be used to recover x{t). However, since 
there are m such functions, this results in more measurements than actually needed. 

We derive the proposed sampling scheme in three steps: First, we consider the problem of compressively 
measuring the vector sequence d[n], whose £th component is given by di[n], where only k out of the m sequences 



de[n\ are non-zero. We show that this can be accompUshed by using the matrix A above and IMV recovery theory. 
In the second step, we obtain the vector sequence d[n] from the given signal x{t) using an appropriate filter bank 
of m analog filters, and sampling their outputs. Finally, we merge the first two steps to arrive at a bank of p < m 
analog filters that can compressively sample x{t) directly. These steps are detailed in the 3 ensuing subsections. 

A. Union of Discrete Sequences 

We begin by treating the problem of sampling and recovering the sequence d[n]. This can be accomplished by 
using the IMV model introduced in Section |IVl Indeed, suppose we measure d[n] with a size p x m matrix A, 
that allows for CS of A;-sparse vectors of length m. Then, for each n, 

y[n]=Ad[n], n e Z. (28) 

The system of ( [281 ) is an IMV model: For every n the vector d[n] is /c-sparse. Furthermore, the infinite set of 
vectors {d[n],n G Z} has a joint sparsity pattern since at most k of the sequences ddn] are non-zero. As we 
described in Section |lVl such a system of equations can be solved by transforming it into an equivalent MMV, 
whose recovery properties are determined by those of A. Since A was designed such that CS techniques will 
work, we are guaranteed that d[re] can be perfectly recovered for each n (or recovered with high probability). The 
reconstruction algorithm is depicted in Fig. [3l Note that in this case the integral in computing Q becomes a sum 
over n: Q = ^^g^yNy^i^^] (we assume here that the sum exists). 

Instead of solving (l28l) we may also consider the Frequency-domain set of equations: 

y(e^'^) = Ad(e-''^), < tj < 27r, (29) 

where y{e^'^),d{e^'^) are the vectors whose components are the frequency responses Y£{e^^) , D i{e^^) . In principle, 
we may apply the CTF block of Fig. [3] to either representations, depending on which choice offers a simpler method 
for determining a basis V for the range of {y(A)}. 

When designing the measurements (l28l ). the only freedom we have is in choosing A. To generalize the class of 
sensing operators we note that d[n] can also be recovered from 

y{e^^) = W{e^'^)Ad{e^^), < w < 27r, (30) 

for any invertible pxp matrix W(e-''^) with elements Wi£{e^^). The measurements of ( [30l l can be obtained directly 
in the time domain as 

p / m \ 

yi[n] = '^Wii[n] * I ^ Afr4[n] I , 1 < i < (31) 

1=1 \r=l / 



where ^i;^^[n] is the inverse transform of Wie{e^'^), and * denotes the convolution operator. To recover d£[n\ from 
y(e-^^), we note that the modified measurements y{e^'^) = W^^(e-''^)y(e-''^) obey an IMV model: 

y(e^'^) = Ad(e^'^), < w < 27r. (32) 

Therefore, the CTF block can be applied to y(e-''^). As in (|3TI) . we may use the CTF in the time domain by noting 
that 

p 

Vi w = X] * ' ^^^^ 

where hit\n\ is the inverse DTFT of Buie?"^), with B(eJ"') = W-^ie?"^). 

The extra freedom offered by choosing an arbitrary invertible matrix W(e-''^) in ( [30l ) will be useful when we 
discuss analog sampling, as different choices lead to different sampling functions. In Section |Vl] we will see an 
example in which a proper selection of W(e-''^) leads to analog sampling functions that are easy to implement. 

B. Biorthogonal Expansion 

The previous section established that given the ability to sample the m sequences di[n\ we can recover them 
exactly from p < m discrete-time sequences obtained via (l30l ) or (OTI ). Reconstruction is performed by applying 
the CTF block to the modified measurements either in the frequency domain (l32l ) or in the time domain (1331 ). The 
drawback is that we do not have access to di[n\ but rather we are given x{t). 

In Fig. |2] and Section Ull we have seen that the sequences d£[n] can be obtained by sampling x{t) with a set 
of functions h£{t) for which M//^(e-''^) of ([T5] ) is stability invertible, and then filtering the sampled sequences 
with a multichannel discrete-time filter 'M.]^\{e^^). Thus, we can first apply this front-end to x{t), which will 
produce the sequence of vectors d[n]. We can then use the results of the previous subsection in order to sense 
these sequences efficiently. The resulting measurement sequences y£ [n] are depicted in Fig. |4l where A is a p x m 
matrix satisfying the requirements of CS in the appropriate dimensions, and W(e-''^) is a size p x p filter bank 
that is invertible a.e. 

Combining the analog filters hi{t) with the discrete-time multichannel filter Mj^^(e-''^), we can express diln] 

as 

de[n] = {ve{t-nT),x{t)), l<£<m,neZ, (34) 

where 

v(w)=M^*^(e^'^^)hH. (35) 
Here v{Lu),h.{uj) are the vectors with £th elements Ve{uj), H£{uj) and (•)"* denotes the conjugate of the inverse. 
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Fig. 4. Analog compressed sampling with arbitrary filters hi{t). 

The inner products in (|34l ) can be obtained by filtering x{t) with the bank of filters vl{—t), and uniformly sampling 
the outputs at times nT. 

To see that ( [34l ) holds, let ci[n] be the samples resulting from filtering x{t) with the m filters v^{t) and uniformly 
sampling their outputs at rate 1/T. From Proposition [U 

c{e^^) = mvA{e^^)d{e^'^). (36) 

Therefore, to establish (l34l ) we need to show that My^(e-''^) = I. Now, from (fTSl ). 

[Mv.(e-)]. = ^ E - - y') 

= i EiM„',(.-)i„ E ^; (y P) 

r=l fcgZ ^ ^ ^ ^ 

= [M^ye^'-)]^[MHA(e^'")]£ = hi, (37) 

where [Q]ir is the irth element of the matrix Q, and [Q]', [Q]i are the ith row and column respectively of Q. 
Therefore, as required, M.yAie^'^) = I. 

The functions {v£{t — nT)} have the property that they are biorthogonal to {a£{t — nT)}, that is 

{vi{t - nT),ai{t - rT)) = bubnr, (38) 

where ba = \ \i I = i, and otherwise. This follows from the fact that in the Fourier domain, (1381 ) is equivalent 



to 



MvAien = I- 



(39) 



Evidently, we can construct a set of biorthogonal functions from any set of functions h£{t) for which WLnAi^'^^) 
is stably invertible, via (l35l ). Note that the biorthogonal vectors in the space A are unique. This follows from the 
fact that if two sets {vj{t)},{vf{t)} satisfy (l38]l, then 



where gi{t) = vl{t) — vf{t). Since {ai{t — nT)} span A, (l40l ) implies that gi{t — mT) lies in A'^ for any i,m. 
However, if both vl{t) and vj{t) are in A, then so is gi{t — mT), from which we conclude that gi{t — mT) = 0. 
Thus, as long as we start with a set of functions hi{t) that span A, the sampling functions Vi{t) resulting from 
(|35]| will be the same. However, their implementation in hardware is different, since hi{t) represents an analog 
filter while M//^(e-^'^) is a discrete-time filter bank. Therefore, different choices of hi{t) lead to distinct analog 
filters. 

C. CS of Analog Signals 

Although the sampling scheme of Fig. |4] results in compressed measurements ye[n\, they are still obtained by an 
analog front-end that operates at the high rate m/T. However, our goal is to reduce the rate at the analog front-end. 
This can be easily accomplished by moving the discrete filters 'MJj\{e^'^), AW(e-''^) back to the analog domain. 
In this way, the compressed measurement sequences yeln] can be obtained directly from x{t), by filtering x{t) 
with p filters S£{t) and uniformly sampling their outputs at times nT, leading to a system with sampling rate p/T. 
An explicit expression for the resulting sampling functions is given in the following theorem. 

Theorem 2: Let the compressed measurements ye[n],l < £ < phe the output of the hybrid filter bank in Fig. |4l 
Then {ye[n]} can be obtained by filtering x{t) of (fTSl l with p filters {s|(— t)} and sampling the outputs at rate 
l/T, where 



Here s{lo), h{Lo) are the vectors with ^th elements Se{Lo), H^{uj) respectively, and the components Ve,{uj) of v(a') = 
mjj\{e^^^)\Y{uj) are Fourier transforms of generators Vi{t) such that {vii{t—nT)} are biorthogonal to {aii{t—nT)} . 



{gi{t - nT),ai{t - rT)) = 0, for all i, i, n, r, 



(40) 



W*(eJ''^^)A*v(a;) 
W*(e^''^^)A*M^*^(e^''^^)h(a;). 



(41) 



In the time domain, 

m p 

= E E E K\-n]A:Mt - nT), (42) 

e=l r=l n£Z 

where Wirl?^] is the inverse transform of [W{e^'^)]ir and 

m 

<t) = E E ^lA-n]h,{t - nT), (43) 

1=1 n& 

where ilJuln] is the inverse transform of \M~^j^{e^^y\ii. 

Proof: Suppose that x{t) is filtered by the p filters Si{t) and then uniformly sampled at nT. From Proposition[T] 
the samples can be expressed in the Fourier domain as 

c(e^''^) = msA{e^'^)d{e^^). (44) 

In order to prove the theorem we need to show that M5A(e-''^) = W(e-''^)A. 
Let 

B(e^'^) = W*(e^'")A*, (45) 

so that s{uj) = B(e^''^^)v(w). Then, 

[yisA[enh = ^ E - y') - y') 



T 

r=l keZ 

= [B*{ennMvA{enh (46) 

where [Q]*, [Q]j are the ith row and column respectively of the matrix Q. The first equality follows from the fact 
that B(e-''^) is 2tt periodic. From (l46l) . 

MsAie^"") = B*(e^''^)MyA(e^''^) = W(e^''^)A, (47) 

where we used the fact that My^(e-''^) = I due to the biorthogonality property. 
Finally, if s{uj) = B(eJ'^^)v(a;), then 

m 

= E E ^^^N^K* - nT), (48) 

e=l n&L 

where 6i^[n] is the inverse DTFT of \Bii{e?^)\ii. Using ( |45] ) together with the fact that the inverse transform of 
Qlnie?^) is (tii\—n\, results in (l42l) . The relation ( |43] ) follows from the same considerations. ■ 



Theorem |2] is the main result which allows for compressive sampling of analog signals. Specifically, starting 
from any matrix A that satisfies the CS requirements of finite vectors, and a set of sampling functions hi{t) for 
which M.HAie-''^) is invertible, we can create a multitude of sampling functions Si{t) to compressively sample the 
underlying analog signal x{t). The sensing is performed by filtering x{t) with the p < m corresponding filters, 
and sampling their outputs at rate 1/T. Reconstruction from the compressed measurements yi[n],l < i < p is 
obtained by applying the CTF block of Fig. [3] in order to recover the sequences di[n]. The original signal x{t) is 
then constructed by modulating appropriate impulse trains and filtering with ai{t), as depicted in Fig. |5] 
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Fig. 5. Compressed sensing of analog signals. The sampling functions Siit) are obtained by combining the blocks in Fig. |4] and are given 
in Theorem |2] 

As a final comment, we note that we may add an invertible diagonal matrix Z(e-^'^) prior to multiplication by 
A. Indeed, in this case the measurements are given by 



y(eJ'^) = W(e-'''^)AZ(e^'^)d(e^'^) = Ad(e- 



(49) 



where d{e^^) has the same sparsity profile as d{e^^). Therefore, d{e^'^) can be recovered using the CTF block. 
In order to reconstruct x{t), we first filter each of the non-zero sequences di[ri\ with the convolutional inverse of 

In this section we discussed the basic elements that allow recovery from compressed analog signals: we first 
use a biorthogonal sampling set in order to access the coefficient sequences, and then employ a conventional CS 
mixing matrix to compress the measurements. Recovery is possible by using an IMV model and applying the CTF 
block of Fig. [3] either in time or in frequency. In practical applications we have the freedom to choose W(e^'^) 
and A so that we end up with analog sampling functions that are easy to implement. Two examples are considered 



in the next section. 



VI. Examples 

A. Periodic Sparsity 

Suppose that we are given a signal x(t) that lies in a SI subspace generated by a{t) so that x{t) = J2n€i, d[n\a{t— 
nT'). The coefficients d[n] have a periodic sparsity pattern: Out of each consecutive group of m coefficients, there 
are at most k non-zero values, in a given pattern. For example, suppose that m = 7,k = 2 and the sparsity profile 
is S = {1, 4}. Then d[n\ can be non-zero only at indices n = l + 7£ or n = 4 + 7£ for some integer £. Decomposing 
d[n\ into blocks d[n] of length m, the sparsity pattern of x{t) implies that the vectors {d[n]} are jointly A;-sparse. 

Since x{t) lies in a SI subspace spanned by a single generator, we can sample it by first prefiltering with the 
filter 

where h{t) is any function such that 4>HA{e^'^) defined by dH) is non-zero a.e. on uj, and then sampling the output 
at rate 1/T', as in Fig. [T] With this choice, the samples c[n\ = {q{t — nT'),x{t)) are equal to the unknown 
coefficients d[n\. We may then use standard CS techniques to compressively sample d[n\. For example, we can 
sample d[n] sequentially by considering blocks d[n] of length m, and using a standard CS matrix A designed to 
sample a /c-sparse vector of length-m. Alternatively, we can exploit the joint sparsity by combining several blocks 
and sampUng them together using MMV techniques, or applying the IMV method of Section |IVl However, these 
approaches still require an analog sampling rate of 1/T' . Thus, the rate reduction is only in discrete time, whereas 
the analog sampling rate remains at the Nyquist rate. Since many of the coefficients are zero, we would like to 
directly reduce the analog rate so as not to acquire the zero values of d[n], rather than acquiring them first, and 
then compressing in discrete time. 

To this end, we note that our problem may be viewed as a special case of the general model ([TSl l with ai{t) = 
a{t — (i — 1)),1 < i < m and di[n] = d[i + nT], 1 < i < m with T = mT'. Therefore, the rate can be reduced 
by using the tools of Section |Vl From Theorem [2j we first need to construct a set of functions {vi{t — nT)} that 
are biorthogonal to {a^(t — nT)}. It is easy to see that 



Vf,{t) =q{t -{£-!)), l<i<m, 



(51) 



with Q{ijj) given by (l50l ) constitute a biorthogonal set. Indeed, with this choice 



T—l 

1 j{i-e)io/T 



T 

r=0 



where we defined 



G(e^'") = ^J2Q*(^- 27rA:) ^ ((^ - 2^A;) . (53) 

From (l50l ). G{e^'^) = 1. Combining this with the relation 

1 

_^g-,2WT^5[^]^ (54) 
r=0 

it follows from §^ that MyA(e^''^) = I. 

We now use Theorem [2] to conclude that any sampling functions of the form s{oj) = W*(e-''^^)A*v(u;) with 

given by (l50b can be used to compressively sample the signal at a rate 
p/T = p/ {mT') < 1/T'. In particular, given a matrix A of size pxm that satisfies the CS requirements, we may 
choose sampling functions 

m 

Si{t)=^A*,q{t -{£-!)), l<i<p. (55) 
e=i 

In this strategy, each sample is equal to a linear combination of several values d[n], in contrast to the high-rate 
method in which each sample is exactly equal d[n\. 

As a special case, suppose that T' = 1, a{t) = 1 on the interval [0, 1] and is zero otherwise. Thus, x{t) is 
piecewise constant over intervals of length 1. Choosing h{t) = a(t) in (l50l ) it follows that q{t) = 1. This is because 
?'aaN = {o-it — n),a{t)) = 6[n\ so that (PaaI^^'^) = 1. Therefore, ve{t) = a{t — {i — 1)) are biorthogonal to 
ae{t). One way to acquire the coefficients d[n\ is to filter the signal x{t) with v{t) = a{t) and sample the output 
at rate 1. This corresponds to integrating x{t) over intervals of length one. Since x(t) has a constant value d[n] 
over the nth interval, the output will indeed be the sequence d[n\. To reduce the rate, we may instead use the 
p < m sampling functions ( |55] ) and sample the output at rate 1/m. This is equivalent to first multiplying x{t) by 
p periodic sequences with period m. Each sequence is piecewise constant over intervals of length 1 with values 
Aj^, 1 < £ < m. The continuous output is then integrated over intervals of length m to produce the samples 
ye[n\. Applying the CTF block to these measurements allows to recover d[n]. Although this special case is rather 



simple, it highlights the main idea. Furthemiore, the same techniques can be used even when the generator a{t) 
has infinite length. 



B. Multiband Sampling 

Consider next the multiband sampling problem [21], [23] in which we have a complex signal that consists of 
at most N frequency bands, each of length no larger than B. In addition, the signal is bandlimited to 27r/T. If 
the band locations are known, then we can recover such a signal from nonuniform samples at an average rate of 
NB /{27r) which is typically much smaller than the Nyquist rate 2tt/T [18], [19], [20]. When the band locations 
are unknown, the problem is much more complicated. In [21] it was shown that the minimum sampling rate for 
such signals is NB/tt. Furthermore, explicit algorithms were developed which achieve this rate. Here we illustrate 
how this problem can be formulated within our framework. 

Dividing the frequency interval [0, 27r/T) into m sections, each of equal length 27r/(mT), it follows that if 
m < 2Tr /{BT) then each band comprising the signal is contained in no more than 2 intervals. Since there are 
bands, this implies that at most 2N sections contain energy. To fit this problem into our general model, let 

Aiioj) = V^A (oj- ^^(^ ~ ^) ] ^ 1 < i < m, (56) 
\ mT J 

where A{oj) is a low-pass filter (LPF) on [0, 27r/(mT)). Thus, Ai{oj) describes the support of the ith interval. 
Since any multiband signal x{t) is supported in the frequency domain over at most 2N sections, x{t) can be 
written as 

m 

X{u:) = Y,M^)Di{u:), (57) 

i=l 

for some Di{uj) supported on the ith interval, where at most 2N functions are nonzero. Since the support of Di{uj) 
has length 2tt /{mT), it can be written as a Fourier series 

Di{uj) = di[n]e-^'^"''^^ = A(e^'"™^). (58) 

Thus, our signal fits the general model ([TSl l. where there are at most 2N sequences di[n] that are nonzero. 

We now use our general results to obtain sampling functions that can be used to sample and recover such signals 
at rates lower than Nyquist. One possibility is to choose hi{t) = ai{t). Since the functions ai{t) are orthonormal 
(as is evident by considering the frequency domain representation), we have that Mj;f^(e-''^) = M^^(e-''^) = I. 
Consequently, the resulting sampling functions are 

m 

Si(t)=5^A:,a,(t). (59) 
1=1 



In the Fourier domain, Si{uj) is bandlimited to 27r/r and piecewise constant with values V mTA*^ over intervals 
of length 2TT/{mT). 

Alternative sampling functions are those used in [21]: 



Si{t) = 6{t-CiT), l<i<p, (60) 

where {cj} are p distinct integer values in the range 1 < Cj < m. Since x{t) is bandlimited, sampling with the 
filters (l60l) is equivalent to using the bandlimited functions 

271- 

S^{Lo) = e-^''^^, 0<u<—. (61) 

To show that these filters can be obtained from our general framework incorporated in Theorem |2j we need to 
choose a. p X m matrix A and an invertible p x p matrix W(e^'^) such that Gi{uj) = Si{uj) where 

g(a;)=W*(e-''-"^^)AV(o;), (62) 

and v(a;) represents a biorthogonal set. In our setting, we can choose Vi{u!) = Ai{u!) due to the orthogonality of 

ai{t). 

Let A be the matrix consisting of the rows a,! < i < p of the m x m Fourier matrix 

A,, = 1 e^2-(^-i)-./™, (63) 



m 



and choose W(e-''^) as a diagonal matrix with ith diagonal element 

Wiie^^) = ^e^'^-'^/'", < cu < 27r. (64) 

From (l62l ). 

m 

G^iu) = Wtie^^^^) AlMu:). (65) 

1=1 

Since Ai,{u)) is equal to \/rnT over the £th interval [{(. - 1)2tt / {raT) 42ti / mT) and otherwise, YJT=i ^lAii'^) 
is piecewise constant with values equal to \/ mTA*^ on intervals of length 27r/(mT). In addition, on the Ah 
interval. 



= {^/VT)e-^^^^^Aie. (66) 
Consequently, on this interval, Gi{uj) is equal to V niTW* {e^'^"^'^)A*^ = e~-'^»'^^. Since this expression does not 



depend on I, Gi{uj) = 5i((j) for < lj < 2tt/T. 

From our general results, in order to recover the original signal x{t) we need to apply the CTF to the modified 
measurements y(e-''^) = W~^(e-''^)y(e-^'^). Since W(e-''^) is diagonal, the DTFT of the ith sequence yi\n\ is given 
by 

Yi{e^'^) = -l=e-^''*'^l'^Yi{e^'^), < a; < 27r. (67) 
yT 

This corresponds to a scaled non-integer delay Cj/m of yi[n]. Such a delay can be realized by first upsampUng 
the sequence yi\n\ by factor of m, low-pass filtering with a LPF with cut-off 7r/m, shifting the resulting sequence 
by Ci, and then down sampling by m. This coincides with the approach suggested in [21] for applying the CTF 
directly in the time domain. Here we see that this processing follows directly from our general framework. 

We have shown that a particular choice of A and W(e^'^) results in the sampling strategy of [21]. Alternative 
selections can lead to a variety of different sampling functions for the same problem. The added value in this 
context is that in [21] there is no discussion on what type of sampUng methods lead to stable recovery. The 
framework we developed in this paper can be applied in this specific setting to suggest more general types of 
stable sampling and recovery strategies. 

VII. Conclusion 

We developed a general framework to treat sampUng of sparse analog signals. We focused on signals in a SI space 
generated by m kernels, where only k out of the m generators are active. The difficulty arises from the fact that 
we do not know in advance which k are chosen. Our approach was based on merging ideas from standard analog 
sampling, with results from the emerging field of CS. The latter focuses on sensing finite-dimensional vectors that 
have a sparsity structure in some transform domain. Although our problem is inherently infinite-dimensional, we 
showed that by using the notion of biorthogonal samphng sets and the recently developed CTF block [38], [21], 
we can convert our problem to a finite-dimensional counterpart that takes on the form of an MMV, a problem 
which has been treated previously in the CS literature. 

In this paper, we focused on sampling using a bank of analog filters. An interesting future direction to pursue 
is to extend these ideas to other sampling architectures that may be easier to implement in hardware. 

As a final note, most of the literature to date on the exciting field of CS has focused on sensing of finite- 
dimensional vectors. On the other hand, traditional sampling theory focuses on infinite-dimensional continuous- 
time signals. It is our hope that this work can serve as a step in the direction of merging these two important areas 
in sampling, leading to a more general notion of compressive sampling. 
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